Libraries
library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
layout(matrix(1:1, nrow=1))
Parameters and
risk
censoredProb <- 0.0002
timeSpan <- 20
timeInterval = 0.01
InitialPopulatoin <- 1000
ContBetaRate_1 <- 0.0002
ContBetaRate_2 <- 0.000001
BinBetaRate_1 <- 0.001
BinBetaRate_2 <- 0.003
BaselineHazard <- ContBetaRate_1
betaRates <- c(BaselineHazard,ContBetaRate_1,ContBetaRate_2,BinBetaRate_1,BinBetaRate_2)
BaseVar <- rep(1,InitialPopulatoin)
ContVar_1 <- runif(InitialPopulatoin)
summary(ContVar_1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0001145 0.2347710 0.4894608 0.4923210 0.7424088 0.9995301
ContVar_2 <- rnorm(InitialPopulatoin,50,10)
ContVar_2[ContVar_2 < 1] <- 1
summary(ContVar_2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 16.41 42.91 49.26 49.62 55.89 83.67
BinVar_1 <- rbinom(InitialPopulatoin,1,0.35)
table(BinVar_1)
## BinVar_1
## 0 1
## 639 361
BinVar_2 <- rbinom(InitialPopulatoin,1,0.35)
table(BinVar_2)
## BinVar_2
## 0 1
## 648 352
dataFeatures <- as.matrix(cbind(BaseVar,ContVar_1,ContVar_2,BinVar_1,BinVar_2))
hazardRate <- as.numeric(dataFeatures %*% betaRates)
summary(hazardRate)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0002355 0.0003632 0.0013097 0.0017651 0.0033525 0.0044680
Getting the events
and time to event
aliveSet <- c(1:InitialPopulatoin)
eventSet <- numeric(InitialPopulatoin)
timetoEvent <- numeric(InitialPopulatoin)
for (time in c(1:(timeSpan/timeInterval)))
{
randProb <- runif(length(aliveSet))
Iscensored <- randProb <= censoredProb
Isevent <- randProb <= (1.0-exp(-hazardRate[aliveSet]))
Iscensored <- Iscensored & !Isevent
eventSet[aliveSet] <- Isevent
timetoEvent[aliveSet] <- time*timeInterval-timeInterval/2
isCensoredOrEvent <- Iscensored | Isevent
aliveSet <- aliveSet[!isCensoredOrEvent]
# cat(length(aliveSet),"(",sum(isCensoredOrEvent),",",sum(Isevent),",",sum(Iscensored),")\n")
}
timetoEvent[aliveSet] <- time*timeInterval + timeInterval/2
summary(timetoEvent)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.005 1.603 5.040 8.283 16.183 20.005
table(eventSet)
## eventSet
## 0 1
## 203 797
pevent <- (1.0-exp(-hazardRate))
summary(pevent)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0002355 0.0003631 0.0013088 0.0017624 0.0033469 0.0044580
simulatedDataFrame <- as.data.frame(cbind(status=eventSet,time=timetoEvent,pevent=pevent,dataFeatures))
RRplots()
plotTimeInterval <- 2.0
hazard <- -log(1.0-simulatedDataFrame$pevent)
hboost <- plotTimeInterval/timeInterval
pvalue <- 1.0-exp(-hboost*hazard)
rdata <- cbind(simulatedDataFrame$status,pvalue)
summary(rdata[,2])
Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04600 0.07007 0.23044 0.26640
0.48855 0.59082
table(simulatedDataFrame$status)
0 1 203 797
RRAnalysisCI <- RRPlot(rdata,atRate=c(0.90,0.80),
timetoEvent=simulatedDataFrame$time,
title="Simulation",
ysurvlim=c(0.00,1.0),
riskTimeInterval=plotTimeInterval)






par(op)
By Risk
Categories
obsexp <- RRAnalysisCI$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
797 |
742.6 |
854.3 |
758.2 |
1.05 |
0.979 |
1.13 |
0.162 |
| low |
190 |
163.9 |
219.0 |
162.4 |
1.17 |
1.010 |
1.35 |
0.034 |
| 90% |
32 |
21.9 |
45.2 |
29.3 |
1.09 |
0.746 |
1.54 |
0.579 |
| 80% |
575 |
529.0 |
624.0 |
567.9 |
1.01 |
0.931 |
1.10 |
0.753 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Cal. Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

Time to Event
Analysis
isevent <- RRAnalysisCI$timetoEventData$eStatus == 1
exptime <- boxplot(RRAnalysisCI$timetoEventData$expectedTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
xlab="Class",
ylab="Time",
main="Expected Time",plot=FALSE)
obstime <- boxplot(RRAnalysisCI$timetoEventData$eTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
xlab="Class",
ylab="Time",
main="Observed Time",
at=exptime$stats[3,],plot=FALSE)
classnames <- attr(RRAnalysisCI$timetoEventData,"ClassNames")
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
Table continues below
| 1Q |
3.26 |
3.89 |
| Median |
8.03 |
8.43 |
| 3Q |
14.18 |
14.26 |
Table continues below
| 1Q |
1.02 |
13.3 |
| Median |
2.35 |
14.8 |
| 3Q |
5.34 |
17.0 |
| 1Q |
11.4 |
1.45 |
| Median |
11.6 |
1.51 |
| 3Q |
11.8 |
3.65 |
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
xlab="Mean Expected Time",
ylab="Mean Observed",
xlim=c(minv,maxv),
ylim=c(minv,maxv),
main="Estimated Time to Event",col=rainbow(length(classnames)))
lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)

Risk
Calibration
op <- par(no.readonly = TRUE)
crdata <- cbind(simulatedDataFrame$status,pvalue,simulatedDataFrame$time)
#calprob <- CalibrationProbPoissonRisk(crdata,timeInterval=plotTimeInterval)
calprob <- CalibrationProbPoissonRisk(crdata)
( 20.72695 , 11020.74 , 1.956533 , 797 , 1588.147 )
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
By Risk
Categories
obsexp <- rrAnalysisTrain$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
797 |
742.6 |
854.3 |
811.7 |
0.982 |
0.915 |
1.05 |
0.623 |
| low |
190 |
163.9 |
219.0 |
173.8 |
1.093 |
0.943 |
1.26 |
0.225 |
| 90% |
32 |
21.9 |
45.2 |
31.4 |
1.019 |
0.697 |
1.44 |
0.858 |
| 80% |
575 |
529.0 |
624.0 |
607.9 |
0.946 |
0.870 |
1.03 |
0.187 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Cal. Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

Time to Event
Analysis
isevent <- RRAnalysisCI$timetoEventData$eStatus == 1
exptime <- boxplot(RRAnalysisCI$timetoEventData$expectedTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
xlab="Class",
ylab="Time",
main="Expected Time",plot=FALSE)
obstime <- boxplot(RRAnalysisCI$timetoEventData$eTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
xlab="Class",
ylab="Time",
main="Observed Time",
at=exptime$stats[3,],plot=FALSE)
classnames <- attr(RRAnalysisCI$timetoEventData,"ClassNames")
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
Table continues below
| 1Q |
3.26 |
3.89 |
| Median |
8.03 |
8.43 |
| 3Q |
14.18 |
14.26 |
Table continues below
| 1Q |
1.02 |
13.3 |
| Median |
2.35 |
14.8 |
| 3Q |
5.34 |
17.0 |
| 1Q |
11.4 |
1.45 |
| Median |
11.6 |
1.51 |
| 3Q |
11.8 |
3.65 |
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
xlab="Mean Expected Time",
ylab="Mean Observed",
xlim=c(minv,maxv),
ylim=c(minv,maxv),
main="Estimated Time to Event",col=rainbow(length(classnames)))
lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)
